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Abstract. The properties of the mean first passage time in a system characterized 
by multiple periodic attractors are studied. Using a transformation from a high 
dimensional space to ID, the problem is reduced to a stochastic process along the 
path from the fixed point attractor to a saddle point located between two neighboring 
attractors. It is found that the time to switch between attractors depends on the 
effective size of the attractors, r, the noise, e, and the potential difference between the 
attractor and an adjacent saddle point as: T — ^exp(^AZ^) ; the ratio between the 
sizes of the two attractors affects /\U. The result is obtained analytically for small 
T and confirmed by numerical simulations. Possible implications that may arise from 
the model and results are discussed. 
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1. Introduction 

The problem of evaluating the mean first passage time (MFPT), has been investigated 
extensively in the context of chemical reactions, dynamical systems etc. More 
recently this problem has been applied to the biological domain, e.g., kinesin walking 
on a microtubule 0], DNA transport through membrane channels [HI, and the time to 
switch between protein folding states 0. 

In this paper the focus is set on dynamical systems characterized by multiple 
periodic attractors. The system may be continuous, however the formal derivation 
is given for discrete time. The problem can be viewed as follows; consider the dynamics 
in a phase space governed by attractors, each represents a dynamic variable. The system 
evolves by changing the amplitudes of these variables. The system being in a particular 
attractor is reflected by the dominance of the associated variable in the current system's 
state. In the absence of noise, the system evolves into one of the available attractors, 
depending on the initial condition. However, when noise is added to the dynamics the 
system may escape from the basin of attraction. The MFPT is then defined as the mean 
time it takes to drive such a system from an attractor to an adjacent saddle point. A 
special case of this problem that arises in the context of a neural network model with 
a feedback loop has been treated previously. This family of model networks is in fact 
a system whose dynamics is influenced by the structure of the attractor space [Zj. A 
geometrically-driven approach to the problem has been proposed in which the MFPT is 
calculated along a ID path embedded in the higher dimensional attractor space. This 
path is a valley connecting the (periodic) fixed-point attractor and the metastable saddle 
point. In periodic systems and the absence of noise, the coupled equations evolve into 
one of the mutually competing fixed points (f.p.'s) in which only the Fourier components 
representing this attractor has a non-vanishing coefficient, i.e. only one of the periodic 
orbits exists asymptotically. The addition of noise generates a perturbation in each of 
the coupled equations. The perturbation can "kick" the system out of the vicinity of 
one stable f.p. so that it escapes to the other f.p. Our motivation is to quantify the 
mean time for such an event to occur. 

In this paper the derivation is extended to arbitrary asymmetric attractors that 
arises from a more general nonlinear dynamical system. An analytical result for the 
MFPT in the limit of weak noise and a weakly non-linear map is derived and confirmed 
numerically. The main reasons for taking these limits are as follows; a high level of 
noise may drive the system too fast from one basin of attraction to another, hence the 
actual dynamics of the system becomes less relevant. The reason for taking weak non- 
linearity is two-fold; one is to facilitate the analytical derivation and allows for certain 
approximations; the other is that most of the characteristics are already revealed in this 
limit. 
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2. The Model 

The family of systems analyzed here can be described by a set of coupled nonlinear 
recurrent equations, generally given by the map: 

x^+i = /.(X„)+a (1) 

where f{X) is a non-linear function, a component of X. The noise term, ^, is a 
Gaussian additive noise distributed according to: 

Without loss of generality and to facilitate the derivations and the graphical 
presentation, we shall focus on the 2D case. In this case the map would be: 

Xn+1 = fl{Xn,yn)+^i 

(3) 

Vn+l = f2{Xn,yn)+^l 

Since the model deals with periodic systems, bounded nonlinear maps are of 
interest. In this case each attractor may be characterized by a dominant Fourier 
component. Generally speaking, the periodic function may be expanded as an infinite 
series, however the demand for weak non-linearity allows for keeping only terms up to 
certain degree, e.g. 3'rd order. A typical such case might be: 

A-^' = g{Y,Al,cos{umt)^ (4) 

where g{) is the nonlinear function and A is the amplitude of a Fourier component. The 
following map (0) is obtained by assuming a bounded nonlinear function that has odd 
terms in the polynomial expansion. Taking only the I'st and 3'rd order terms of g leads 
to: 

Xn+1 = (1 + Ti)xn [1 " ax^ - byl] 

(5) 

Vn+l = (1 + T'2)yn [1 " ay^ - bxl] 

where < a, 6 < 1 are constants and |xo| < 1, \yo\ < 1. 

The key point is our ability to identify a low dimensional discrete dynamics that 
describes the evolution of the system, i.e., the amplitude of the dominant Fourier 
components in the asymptotic meta-stable state (w/o noise). 

Using the following variables rescaling: 



X = s/T^x y = y/ny 



(6) 
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gives rise to the potential function 



2 A Ti T2 2 



(7) 



where the coupled equations are obtained via 



Vn+l = Vn 



T2K . 



Keeping only terms with the small parameter r to first order, the following transformed 
equations are obtained: 



Xn+l 



Xn Tl Xfi '-\~ CL -|- bXfiyfi 

Vn+l = Vn - T2 [^Vn + a^vl + ^^n^n 



(8) 



The fixed points of the dynamics (assuming Tj <^ 1) are given (for the rescaled 
variables) by 



aT2 











aTi 



± 



h^T2 



± 



,7-2 
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where the last (four) f.p. are in fact saddle points and the first (four) are stable f.p. An 
example of the phase portrait in this case is given in figure ^ 

To visualize the potential surface an example is depicted in figure |2| for clarity 
the negative potential is shown - more 'red' means more negative potential (e.g., the 
zero F.P. is not stable). The variables shown are rescaled and the actual values of the 
potential are not important. It is apparent that there exists a path from a non-zero F.P. 
to adjacent S.P. via a valley (here a hill). 

As mentioned above, an important step in the derivation is the conjecture that 
the most probable trajectory from a f.p. to one of its nearest s.p. defines the properties 
of the MFPT. In the next section the properties of the rescaled ID noisy map with a 
potential given by equation dare derived. By fixing the initial and final points the path 
is unique, hence, the projection of near by trajectories on the path are treated and the 
noise term is rescaled appropriately. 



3. Escape from a meta-stable attractor 

To illustrate the idea behind the following derivation let us assume a scenario where 
the initial condition is ?/ = , x = x*, i.e. one of the f.p.'s of the dynamics (refer 
to figure Q). Since the line connecting this f.p. and the saddle point is a valley (in 
the potential space), it may be assumed that the most probable escape route is along 
this line (or its mirror through the x-axis, i.e. the line connecting the f.p. with the 
saddle point (xj^, ). This argument can be understood by rotating each noise term 
tangent and perpendicular to the path. The perpendicular term decays fast due to the 
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restoring force; hence, the crucial step that the dynamics is mainly ID is conjectured. 
It should be noted that this conjecture is not limited to our 2D example; rather, a set 
of attractors embedded in a higher dimensional space should exhibit the same property, 
namely, that the most probable escape route from the domain of attraction is through 
the valley connecting this attractor and a saddle-point residing between two adjacent 
f.p.'s. Therefore, with the assumption of weak noise and r <C 1 the map can be reduced 
into one dimension, on that path; hence, a ID noisy map is obtained: 



Sn+l = Sn- rU'iSn) + , (10) 

where s defines the path. The noise term now is the tangential projection of the noise 
vector on the path. The explicit form of the path is not crucial as long as smoothness and 
monotonicity are guaranteed; since it is assumed that there are no extremums between 
the f.p. and the s.p. this assumption holds. 

This type of a ID equation has been investigated by several authors |2llHlin! for the 
case of small non-linearity, namely, the class of map functions with the property that 
f{s) deviates only weakly from the identity map: 



/(^) = ^-r^ , r«l . (11) 

in analogy with our ID map (equation [TUI) . 

Assume that the process described in equation [TUI is defined in (— cx), oo) and let us 
define the random variable t{s), the first passage time from the interval I = [SP~ , SP~^] 
,by: 

i=mm{n : > s+} , (12) 

i.e. the first time the process hit one of the boundaries, where SP^ are the saddle points 
defined above, and is the value of s at the saddle point. The MFPT, t{s), starting 
from a point in / is given by: 

t{s) =< i{s) >= E[i\ So = s] . (13) 
It was shown (e.g., jH]) that the MFPT can be written as 



t{s) - 1 = y Piz\s)t{z)dz , (14) 

where P(z\s) denote the transition probability to go from s„ = s to = z in a 
single step. 

In the next subsections the probability density function and the MFPT are derived. 
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3.1. Probability density function 

Assume the noise terms, ^„ , are Gaussian distributed and mutual independent, hence, 
the stochastic process s„ is a Markov process. The distribution is given by equation El 
with a rescaled amphtude e. The transition probabihty from s„ = 2; to = s in one 
step is 

P{s\z) = p{s-f{z)) , (15) 

with p(^) given by equation |21 The probabihty density, Qn{s), to find the system in s 
after n steps evolves according to 



/oo 
P{s\z)Qn{z)dz . 
-00 



(16) 



In a similar way, one can define the conditional probability Pn{s\z) to get from z to 
[s, s + ds] in n steps. The invariant probability density Q{s) is the solution of equation 
Uni Plugging eq. HSl the invariant density obey 

Q{s) = J P{s\z)Q{z)dz = {27re)-^/^ J exp i^-^Js- z + Tie {z)]'^^Q{z)dz .(17) 

This integral can not be solved in general. In the case of weak non-linearity, 
e.g. our map with r ^ 1, one can solve the integral approximately by setting 
u = (2e)^^/^ (s — 2; + rW(2)) and neglecting C(t^) terms in the expansion of z{u). 
Alternatively, one can assume a WKB-type solution (e.g., jTU]) to solve equation [T7I for 
weak noise, e <^ 1, 



Q(s) = A^(s)exp 



where $(s) is a potential function and N{s) is a normalization factor. Inserting this 
type of solution back in equation [T71 one obtains: 



N{s) = (27re)~i/2 / N{z)exp 



-^W'{.s,z) 



dz 



(19) 



where 



W\s,z)=2mz)-^s)] + [s-fiz)f , (20) 

whereas in our model f{z) = z — tU'{z). The demand for (semi) positivity of 
implies that the potential $(s) is a Lyapunov function of the noise- free dynamics. This 
is evident by setting s = f{z) in equation 1201 To obtain an equation for the potential 
$(s) one can derive a transformation from the old variable, z, to the new one, W. It is 
sufficient and necessary that dW{s, z)/dz 7^ to write 



N{s) = (27re)-i/2 / N{z{W,s)) [dW{s, z)/dzl=,^w,s)] ' 



exp 



1 . 
— W 
2e 



dW ,(21) 
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where z(W. s) is a solution of equation QUI Assuming N{z{y, s)) [dW/dz\z=z{w,s)] ^ is 
smooth in a neighborhood of (9(e^/^) around W = 0, the integral can be approximated 
for weak noise by 



N{s)=N{z{y = 0,s))[dW/dz\z=,^w=o,s)] • (22) 

Assuming the variable transformation remains valid ( dW/dz ^ ) for the case 
dW^{s, z)/dz = , must vanish as well. Therefore, differentiating QUI leads to an 
equation for s. Plugging = and the equation for s back to equation 1201 the desired 
relation for $(s) is obtained: 

^'['l . +^-rU'iz)]+U ^ ^'['lX = , (23) 



1-tU"{z) ' 7 2 \1-tU"{z_ 

where the map / has been written explicitly and lA"{z) is a second order derivative w.r.t. 
z. The solution of this equation for small non-linearity, r <^ 1, gives finally 

g(,) = iVexp^^^ , (24) 
e 

where the pref actor is constant up to corrections of O(t^). 



3.2. The MFPT 

Turning now to the analysis of the MFPT, consider the random variable t{s) - the 
random time to leave the domain of attraction, /, for the first time, starting from a 
point s E I. The probabihty that after n steps the process does not leave the domain / 
is given by 

Q{n\s) = jPn{z\s)dz , (25) 

where Pn{z\s) is the probability to get from s to z in exactly n steps, obeying a similar 
equation asUHl 

Pn+Ms) = jpMy)PMs)dy , (26) 

where this recursive relation uses the one-step probability Pi. Since Q{n\s) is a 
decreasing function of the discrete time n, the probability that an exit occurs exactly 
after n steps is simply Q{n\s) — Q{n + l|s) , hence, the first moment of our random 
variable, i{s), is 

oo oo 

t{s) = {i{s)) = Y,in + l){Q{n\s)-Q{n + l\s)) = J2Qin\s) . (27) 

n=0 n=0 
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Summing equation over n and taking the integral inside the domain / only, one 
obtains: 

oo „ oo 

^g(n + = / ^Q(n|s)Pi(2|s)dz . (28) 

n=0 n=0 

On noting that the l.h.s. equals t(s) — 1 ( Q{0\s) = 1 V s G / ) this equation becomes 
equation [T31 Multiplying equation HH by the invariant probability density Q{s) and 
integrating over the domain as described by equation O i.e., I = [SP~ , SP'^] , it is 
found that: 

/ Q{s)ds=\ / ds+ ds]Q{s) / Pi{z\s)t{z)dz . (29) 

JsP- JsP+ J JsP- 

Under the assumption of weak noise e <^ 1, the function t{s) is nearly constant 
inside the domain of attraction. Fluctuations occur mainly near the boundary. The 
reason is that only close to the boundary may one have a finite probability to jump over 
the boundary in a small number of steps. Therefore, it was suggested (see 131^]) that 
this function can be written as a product of a constant value, and a boundary layer 
function: 



t{s) = Th{s) , h{s*) = 1 , (30) 

where s* is the f.p. The boundary layer extends a distance of order e^/^ around s = s^, 
and the scaled boundary layer function h{s) may be written as: h{s) = h({2eY^'^s). 

Using the boundary layer function h (equation IHUj) to account for the fact that for 
weak noise t(s) is almost constant inside the domain, except near the boundaries, MFPT 
(starting from the f.p.) is given by: 

(l^L ds + J^p+dsj Q{s) Jgpl Pi{z\s)h{z)dz 
Jsp- Q{s)ds 

where the first integral in the numerator is the inverse MFPT from the negative 
exit and the second from the positive. Since the saddle points are symmetric w.r.t. 
the f.p. the integrals equal. The invariant density is peaked around the stable f.p., 
hence, the steepest descent approximation may be used for both the numerator and the 
denominator of equation |^ Expanding the map near the f.p., 

f{z* + Az) ^z' + Azf'iz') = z* + Az{l - tU"{z*)) , (32) 

where z* denotes the f.p., and using the Gaussian approximation of equation [TTl the 
denominator of equation |^ becomes 



SP+ 

Qis)ds ^ Q{FP) 

sp- 



ire 



tU"{FP) 



1/2 



(33) 



Mean First Passage Time in Periodic Attractors 9 

where the weak non-hnearity assumption, r ^ 1, is used. The numerator is evaluated 
in a similar manner yielding 

/ ds+ ds]Q{s) Pi{z\s)h{z)dz^ e^/^G{T)Q{SP) ,(34) 

J-oc, JsP+ J JsP- 

With G(r) = (^)V2 + 0(r) . 

Combining the terms, the main result is obtained: 

T = ^^^^=^e^p^{U{SP)-U{FP)) , (35) 

where C is a constant. Note that the time T is rescaled w.r.t. a cycle of the attractor. 

Plugging the general form of the f.p.'s of the dynamics (equation El) in equation 
[71 an explicit expression for the potential terms in equation EHl as a function of the 
variables of the model, ti,T2, a, b, is obtained as follows: 

^(SP) = - "f ^"r'tl" KiFP) = (--^) (36) 

where t\2 = ti/t2\ the parentheses on the r.h.s. are given for the second f.p. Note that 
Ti,T2 are given now w.r.t. a common variable r, e.g., T2 = t and ri = T12T. The special 
case Ti = T2 = T reduces to (U{SP) — U{FP)) = {b — a) / (4a(a + 6)), and for the choice 
a = 1/4, 6 = 1/2 used below for simulations, AU = 1/3. 



3.3. Numerical simulations 

To confirm the theoretical results extensive numerical simulations of the model were 
performed. Equations |H1 were used to evaluate the statistical properties of the MFPT. 
For each trial, a line that passes through the s.p. and perpendicular to the line connecting 
the f.p. and the s.p. was constructed. The time to hit this line starting from the f.p. was 
collected in 500 — 1000 trials for each choice of the parameters. The model parameters 
used in the figure were a = 1/4,6 = 1/2 (although other values were tested as well). 
Typically the attractor and noise parameters were chosen such that r/e G (2. .7), and 
0.8 < ri2 < 1.2. 

To demonstrate the scaling properties of the results, equations EH- ESI the logarithm 
of the MFPT is evaluated and the following quantity is depicted in the next figure: 

<lnT>+lnr oc r/e 

The slope of each series of simulations corresponds to a specific choice of ri2,T, e. The 
simulation results shown in figure El are in a good agreement with the predicted values 
calculated using equation EHl see table [H 
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Ti2 


2(W(SP) - U(FP)) 


Simulation 


0.8 


0.3 


0.32 ± 0.02 


0.9 


0.47 


0.45 ± 0.03 


1.0 


0.66 


0.64 ± 0.03 


1.1 


0.88 


0.90 ± 0.03 


1.2 


1.09 


1.11 ± 0.03 



Table 1. Predicted vs. simulated potential difference values 
4. Discussion 

The theory for the mean first passage time to escape a periodic attractor defined by 
the particular fixed point of the dynamic equations was developed. To facilitate the 
analytical investigation, proper variables transformation were identified that decompose 
the potential function of the system. The reduced variables are closely related to the 
amplitude of the solution. The noiseless system relaxes to one of the stable non-zero 
solution (above bifurcation). Adding noise to the dynamics perturbs this solution and 
enables possible escape. 

One of the key points in the solution is the conjecture that the most probable 
escape route is essentially one dimensional along the path from the attractor to one of 
the adjacent saddle points. This enabled us to reduce the dimensionality of the dynamics 
into a ID process. Taking the limit of small noise and not too far from the bifurcation 
the theory developed for noise driven discrete dynamical systems has been applied. The 
results resemble those obtained in systems with potential barrier undergoing a tunneling 
in the sense that the escape time has a polynomial prefactor and a leading exponential 
term. 

Simulations of the system with two variables have shown that the theory developed, 
and especially the reduction to a ID flow, are applicable and provide a good prediction 
of the exponential term, as well as the polynomial prefactor. 

It should be noted that the results are applicable to non periodic systems as well, 
as long as one can transform the dynamic equations to a form similar to equations 
i.e., up to third order. Also, systems with higher order terms may exhibit a very 
similar behavior assuming these terms do not dominate the solutions and the discussion 
is restricted to solutions not too far from bifurcation, i.e. r ^ 1. 

The implications that arise from this results may have practical application in fields 
such as optimization, minimization problems, numerical analysis, parameter estimation 
in high dimensional complex systems (e.g., neural networks) and more. The common 
ground to these problems is the existence of (usually) high dimensional attractor space 
in which some dynamics is imposed, e.g., obtaining a minimal energy solution to 
optimization problem via some iterative algorithm. This algorithm evolves the solution 
preferably to a global minimum, however it may get trapped in local minima. To escape 
these attractors, one occasionally perturbs the solution to escape from the solution etc. 
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It is of practical importance to assess the time to escape (affect simulation time) and 
possibly the required amplitude and duration of perturbation. All these issues are left 
for future research. In this context it is interesting to mention a work by Baronchelli and 
Loreto ^T] on mean first passage time in graphs, where the focus is on the complexity 
of evaluating the MFPT on a particular node of a graph. 
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Figure 1. Phase portrait of a 2D map witii n = 0.02, T2 = 0.025, a = 0.2, b = 0.4. 
'SP' denotes a saddle point and 'FP' a fixed point. Arrows shows the direction of the 
flow. 




Figure 2. Negative potential surface of tlie function given in equation[7| The insert is 
the contour plot. The parameters used were ^ = 3/4, a = 0.5, b = 0.25. 'SP' denotes 
a saddle point and 'FP' a fixed point. 
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Figure 3. Scaling of the average logarithm of the escape time in the 2D model. The 
dashed lines are the linear regression of the corresponding data points. The standard 
deviation of each point is roughly the same as the symbol, hence omitted for clarity. 
The legend shows for each ri2 — ti/t2 the estimated slope in parentheses, see also 
table □ 



